// gcc H2OToFile.c -o H2O.exe -lglut -lm -lGLU

#include <GL/glut.h>
#include <math.h>
#include <stdio.h>

#include <GL/glut.h>
#include <math.h>
#include <stdio.h>

#define GR 2.0
#define H  1.0
#define G  0.1

#define N 64
#define M 4
#define CUBE_ROOT_N 4

#define SEPERATION 3.0
#define POSITION_START_VALUE (-(3.0*SEPERATION*CUBE_ROOT_N/2.0 - 3.0*SEPERATION/2.0))

//#define SEPERATION 2.0
//#define POSITION_START_VALUE -2.0
 
#define SCALE (8.0*(1.5*CUBE_ROOT_N - 1.5))
#define XMax SCALE
#define XMin -SCALE

#define YMax SCALE
#define YMin -SCALE

#define ZMax SCALE
#define ZMin -SCALE

#define XWindowSize 1000
#define YWindowSize 1000

//Frustum settings
#define LEFT  -0.2
#define RIGHT  0.2
#define UP     0.2
#define DOWN  -0.2
#define NEAR   0.2
#define FAR   10.0

//LookAT settings
#define EYE_AT_X   2.0
#define EYE_AT_Y   2.0
#define EYE_AT_Z   2.0
#define LOOK_AT_X  0.0
#define LOOK_AT_Y  0.0
#define LOOK_AT_Z  0.0
#define UP_X       0.0
#define UP_Y       1.0
#define UP_Z       0.0

#define PI 3.141592654

#define DRAW 1000

#define DT 0.0002
#define TSTART 0.0
#define TSTOP 400.0

//Stuff to draw spheres
#define c 0.525731
#define d 0.850651

//Globals
float NormalizedVerticies[12][3] = {{-c , 0.0,   d}, {c, 0.0, d}, {-c, 0.0, -d}, {c, 0.0, -d},
			  	    {0.0, d  ,   c}, {0.0, d, -c}, {0.0, -d, c}, {0.0, -d, -c},
			            {d  , c  , 0.0}, {-d, c, 0.0}, {d, -c, 0.0}, {-d, -c, 0.0}};
			            
float NormalVectors1[20][3], NormalVectors2[20][3], NormalVectors3[20][3];

static GLuint tindices[20][3] = {{1,4,0} , {4,9,0} , {4,5,9} , {8,5,4} , {1,8,4},
				 {1,10,8}, {10,3,8}, {8,3,5} , {3,2,5} , {3,7,2},
				 {3,10,7}, {10,6,7}, {6,11,7}, {6,0,11}, {6,1,0},
				 {10,1,6}, {11,0,9}, {2,11,9}, {5,2,9} , {11,2,7}};
				 
float g_total_mass;


float WorldToScreen(float coordinate, float min, float max)
{
	return( -1.0 + 2.0*(coordinate - min)/(max - min));
}

//sets vectors that are normal to each of the 20 faces that makeup a sphere
void InitNormalVectors()
{
	int j, k;
	for(j = 0; j < 20; j++)
	{
		for(k = 0; k < 3; k++)
		{
			NormalVectors1[j][k] = NormalizedVerticies[tindices[j][0]][k];
			NormalVectors2[j][k] = NormalizedVerticies[tindices[j][1]][k];
			NormalVectors3[j][k] = NormalizedVerticies[tindices[j][2]][k];
		}
	}
}

void DrawLowResSphere(float x, float y, float z, float Radius, float Red, float Green, float Blue)
{
	int i;
	float VectorData1[20][3], VectorData2[20][3], VectorData3[20][3];

	for(i = 0; i < 20; i++)
	{
		VectorData1[i][0] = WorldToScreen(Radius*NormalVectors1[i][0] + x, XMin, XMax);
		VectorData1[i][1] = WorldToScreen(Radius*NormalVectors1[i][1] + y, YMin, YMax);
		VectorData1[i][2] = WorldToScreen(Radius*NormalVectors1[i][2] + z, ZMin, ZMax);
		VectorData2[i][0] = WorldToScreen(Radius*NormalVectors2[i][0] + x, XMin, XMax);
		VectorData2[i][1] = WorldToScreen(Radius*NormalVectors2[i][1] + y, YMin, YMax);
		VectorData2[i][2] = WorldToScreen(Radius*NormalVectors2[i][2] + z, ZMin, ZMax);
		VectorData3[i][0] = WorldToScreen(Radius*NormalVectors3[i][0] + x, XMin, XMax);
		VectorData3[i][1] = WorldToScreen(Radius*NormalVectors3[i][1] + y, YMin, YMax);
		VectorData3[i][2] = WorldToScreen(Radius*NormalVectors3[i][2] + z, ZMin, ZMax);
	}
	glColor3f(Red,Green,Blue);
	glBegin(GL_TRIANGLES);
	for(i = 0; i < 20; i++)
	{
		glNormal3fv(NormalVectors1[i]);
		glVertex3fv(VectorData1[i]);
		glNormal3fv(NormalVectors2[i]);
		glVertex3fv(VectorData2[i]);
		glNormal3fv(NormalVectors3[i]);
		glVertex3fv(VectorData3[i]);
	}
	glEnd();
}

void set_initail_conditions(float *pos, float *vel, float *mass)
{
	int i,j,k,l;
	float startx, starty, startz;
	
	startx = 0.0 - (float)(N*SEPERATION)/2.0;
	starty = 0.0 - (float)(N*SEPERATION)/2.0;
	startz = 0.0 - (float)(N*SEPERATION)/2.0;

	l = 0;
	for(i = 0; i < 3*CUBE_ROOT_N; i += 3)
	{
		for(j = 0; j < 3*CUBE_ROOT_N; j += 3)
		{
			for(k = 0; k < 3*CUBE_ROOT_N; k += 3)
			{
				pos[l    ] = POSITION_START_VALUE + i*SEPERATION;
				pos[l + 1] = POSITION_START_VALUE + j*SEPERATION;
				pos[l + 2] = POSITION_START_VALUE + k*SEPERATION;
				l += 3;
			}
		}
	}
	
	for(i = 0; i < 3*N; i++)
	{
		vel[i] = 0.0;		
	}
	
	for(i = 0; i < M; i++)
	{
		mass[i] = 1.0;		
	}
	
	for(i = M; i < N; i++)
	{
		mass[i] = 1.0;		
	}
	
	g_total_mass = 0.0;
	for(i = 0; i < N; i++)
	{
		g_total_mass += mass[i];		
	}
}

void draw_picture(float *pos)
{
	int i;

	glClear(GL_COLOR_BUFFER_BIT);
	glClear(GL_DEPTH_BUFFER_BIT);
	glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
	
	for(i = 0; i < 3*N; i += 3)
	{
		DrawLowResSphere(pos[i], pos[i + 1], pos[i + 2], 0.75, 0.0, 0.0, 1.0);
	}
	glutSwapBuffers();
	glFlush();
}

void reset_center_of_mass(float *pos, float *mass)
{
	int i,j;
	float center_of_mass_x, center_of_mass_y, center_of_mass_z;

	center_of_mass_x = 0.0;
	center_of_mass_y = 0.0;
	center_of_mass_z = 0.0;
	
	j = 0;
	for(i = 0; i < 3*N; i += 3)
	{
		center_of_mass_x += pos[i]*mass[j];
		center_of_mass_y += pos[i + 1]*mass[j];
		center_of_mass_z += pos[i + 2]*mass[j];
		j++;
	}
	center_of_mass_x /= g_total_mass;
	center_of_mass_y /= g_total_mass;
	center_of_mass_z /= g_total_mass;

	for(i = 0; i < 3*N; i += 3)
	{
		pos[i    ] -= center_of_mass_x;
		pos[i + 1] -= center_of_mass_y;
		pos[i + 2] -= center_of_mass_z;
	}
}

__global__ void get_forces(float *pos, float *acc, float *mass)
{	
	int m  = blockDim.x;
	
	int ix = 3*blockIdx.x;
	int iy = 3*blockIdx.x + 1;
	int iz = 3*blockIdx.x + 2;
	
	int jx = 3*threadIdx.x;
	int jy = 3*threadIdx.x + 1;
	int jz = 3*threadIdx.x + 2;
	
	float dx, dy, dz, d2, inv_d2, force_mag;
	
	__shared__ float sh_pos[3*N];
	__shared__ float sh_force[3*N];
	__shared__ float sh_mass[N];
	
	if(threadIdx.x < N)
	{
		sh_pos[jx] = pos[jx];
		sh_pos[jy] = pos[jy];
		sh_pos[jz] = pos[jz];
		sh_mass[threadIdx.x] = mass[threadIdx.x];
	}
	__syncthreads();
	
	if(threadIdx.x < N && threadIdx.x != blockIdx.x)
	{
		dx = sh_pos[jx] - sh_pos[ix];
		dy = sh_pos[jy] - sh_pos[iy];
		dz = sh_pos[jz] - sh_pos[iz];
	
		d2 = (dx*dx + dy*dy + dz*dz);
		inv_d2 = 1.0/d2;
		force_mag = -G/(d2*d2*d2) + H/(d2*d2) + GR*sh_mass[blockIdx.x]*sh_mass[threadIdx.x]/d2;
		//force_mag = inv_d2*(sh_mass[blockIdx.x]*sh_mass[threadIdx.x]*GR + inv_d2*(H - inv_d2*G))/sqrt(d2);
		sh_force[jx] = force_mag*dx;
		sh_force[jy] = force_mag*dy;
		sh_force[jz] = force_mag*dz;
	}
	if(threadIdx.x == blockIdx.x)
	{
		sh_force[ix] = 0.0;
		sh_force[iy] = 0.0;
		sh_force[iz] = 0.0;
	}
	__syncthreads();
	
	while(m > 1)
	{	
		m = m/2;
		if(threadIdx.x  < N)
		{
			sh_force[jx] += sh_force[jx + m];
			sh_force[jy] += sh_force[jy + m];
			sh_force[jz] += sh_force[jz + m];
		}
		__syncthreads();
	}	
	__syncthreads();
	
	acc[ix] = sh_force[0]/sh_mass[blockIdx.x];
	acc[iy] = sh_force[0]/sh_mass[blockIdx.x];
	acc[iz] = sh_force[0]/sh_mass[blockIdx.x];
}

void nBody()
{
	int i;
	int tdraw;
	float time;
	float *pos_CPU, *vel_CPU, *acc_CPU, *mass_CPU; //CPU pointers
	float *pos_GPU, *vel_GPU, *acc_GPU, *mass_GPU; //GPU pointers
	dim3 dimBlock, dimGrid;
	
	//Allocate Device (GPU) Memory, & allocates the value of the specific pointer/array
	cudaMalloc(&pos_GPU,3*N*sizeof(float));
	cudaMalloc(&vel_GPU,3*N*sizeof(float));
	cudaMalloc(&acc_GPU,3*N*sizeof(float));
	cudaMalloc(&mass_GPU,N*sizeof(float));

	//Allocate Host (CPU) Memory
	pos_CPU = (float*)malloc(3*N*sizeof(float)); //(float*) to prevent from being a void
	vel_CPU = (float*)malloc(3*N*sizeof(float));
	acc_CPU = (float*)malloc(3*N*sizeof(float));
	mass_CPU = (float*)malloc(N*sizeof(float));

	set_initail_conditions(pos_CPU, vel_CPU, mass_CPU);

	draw_picture(pos_CPU);
	
	//Threads in a block
	dimBlock.x = N;
	dimBlock.y = 1;
	dimBlock.z = 1;
	
	//Blocks in a grid
	//dimGrid.x = (N - 1)/dimBlock.x + 1; // makes enough blocks to finish task
	dimGrid.x = N;
	dimGrid.y = 1;
	dimGrid.z = 1;

	tdraw = 0;
	time = TSTART;

	while(time < TSTOP)
	{
		cudaMemcpy(pos_GPU, pos_CPU, 3*N*sizeof(float),cudaMemcpyHostToDevice);
		cudaMemcpy(mass_GPU, mass_CPU, N*sizeof(float),cudaMemcpyHostToDevice);

		get_forces<<<dimGrid, dimBlock>>>(pos_GPU, acc_GPU, mass_GPU);

		cudaMemcpy(acc_CPU, acc_GPU, 3*N*sizeof(float),cudaMemcpyDeviceToHost);

		//Updates of velocity and position
		for(i = 0; i < 3*N; i += 3)
		{
			vel_CPU[i]     += DT*acc_CPU[i]     ;
			vel_CPU[i + 1] += DT*acc_CPU[i + 1] ;
			vel_CPU[i + 2] += DT*acc_CPU[i + 2] ;

			pos_CPU[i]     += DT*vel_CPU[i];
			pos_CPU[i + 1] += DT*vel_CPU[i + 1];
			pos_CPU[i + 2] += DT*vel_CPU[i + 2];
		}

		time += DT;
		tdraw++;

		if(tdraw == DRAW)
		{
			reset_center_of_mass(pos_CPU, mass_CPU);
			draw_picture(pos_CPU);
			tdraw = 0;
		}
	}
}

void Display(void)
{
	InitNormalVectors();

	gluLookAt(EYE_AT_X,EYE_AT_Y, EYE_AT_Z, LOOK_AT_X, LOOK_AT_Y, LOOK_AT_Z, UP_X, UP_Y, UP_Z);

	nBody();
}

void reshape(int w, int h)
{
	glViewport(0, 0, (GLsizei) w, (GLsizei) h);

	glMatrixMode(GL_PROJECTION);

	glLoadIdentity();

	glFrustum(LEFT, RIGHT, DOWN, UP, NEAR, FAR);

	glMatrixMode(GL_MODELVIEW);
}

int main(int argc, char** argv)
{
	glutInit(&argc,argv);
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_DEPTH | GLUT_RGB);
	glutInitWindowSize(XWindowSize,YWindowSize);
	glutInitWindowPosition(0,0);
	glutCreateWindow("H2O");
	GLfloat light_position[] = {1.0, 1.0, 1.0, 0.0};
	GLfloat light_ambient[]  = {0.0, 0.0, 0.0, 1.0};
	GLfloat light_diffuse[]  = {1.0, 1.0, 1.0, 1.0};
	GLfloat light_specular[] = {1.0, 1.0, 1.0, 1.0};
	GLfloat lmodel_ambient[] = {0.2, 0.2, 0.2, 1.0};
	GLfloat mat_specular[]   = {1.0, 1.0, 1.0, 1.0};
	GLfloat mat_shininess[]  = {10.0};
	glClearColor(0.0, 0.0, 0.0, 0.0);
	glShadeModel(GL_SMOOTH);
	glColorMaterial(GL_FRONT, GL_AMBIENT_AND_DIFFUSE);
	glLightfv(GL_LIGHT0, GL_POSITION, light_position);
	glLightfv(GL_LIGHT0, GL_AMBIENT, light_ambient);
	glLightfv(GL_LIGHT0, GL_DIFFUSE, light_diffuse);
	glLightfv(GL_LIGHT0, GL_SPECULAR, light_specular);
	glLightModelfv(GL_LIGHT_MODEL_AMBIENT, lmodel_ambient);
	glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
	glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess);
	glEnable(GL_LIGHTING);
	glEnable(GL_LIGHT0);
	glEnable(GL_COLOR_MATERIAL);
	glEnable(GL_DEPTH_TEST);
	//glutKeyboardFunc(KeyPressed);
	//glutKeyboardUpFunc(KeyUp);
	//glutMouseFunc(Mouse);
	glutDisplayFunc(Display);
	glutReshapeFunc(reshape);
	//DoWork();
	//glutIdleFunc(IdleFunction);
	//SetUpMenus();
	glutMainLoop();
	return 0;
}






